%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pyBSS as pb
import scipy.linalg as lng
import pyfits as pyf
import copy as cp
import BSS_utilities as bu
from imp import reload
from redwave import RedWave
#
fsize = 12
vcol = ['mediumseagreen','crimson','steelblue','darkmagenta','burlywood','khaki','lightblue','darkseagreen','deepskyblue','forestgreen','gold','indianred','midnightblue','olive','orangered','orchid','red','steelblue']
font = {'family' : 'normal',
'weight' : 'bold',
'size' : fsize}
plt.rc('font', **font)
n = 2 # number of sources
m = 2 # number of observations
t = 1024 # number of samples per observation
# The sources will be :
# - Gaussian is SType = 1
# - Uniform is SType = 2
# - Sparsely-distributed is SType = 3
SType = 3
# noise level in dB, it is basically negligible in this case
noise_level = 120
X,A,S = pb.MakeMixture(SType=SType,n=n,m=m,t=t,noise_level = noise_level)
plt.plot(X.T) # plotting the mixtures
plt.title("Two mixtures of sources")
fig = plt.gcf()
fig.set_size_inches(12, 8)
I.1 - Looking at scatter plots¶
A key tool to better understand how BSS methods work is the scatter plot of the data or sources. Basically, it amounts to representing samples of a given signal (e.g. a source) as a function of another. From the viewpoint of statistics, this is an estimate of the joint probability density function of these signals. It looks like this:
vrange = 0.5*np.max(np.sqrt(np.sum(X*X,axis=0)))
plt.subplot(121)
plt.plot(X[0,:],X[1,:],'.',color='crimson')
plt.xlabel("Data $x_1$")
plt.ylabel("Data $x_2$")
G = A.T
plt.arrow( 0, 0,1.2*vrange*G[0,0], 1.2*vrange*G[0,1], fc="k", ec="k",linewidth=2,head_width=0.1*vrange, head_length=0.1*vrange)
plt.arrow( 0, 0,1.2*vrange*G[1,0], 1.2*vrange*G[1,1], fc="k", ec="k",linewidth=2,head_width=0.1*vrange, head_length=0.1*vrange)
plt.title("Scatter plot of the data")
plt.axis([-2.5*vrange, 2.5*vrange, -2.5*vrange, 2.5*vrange])
plt.subplot(122)
plt.plot(S[0,:],S[1,:],'.',color='steelblue')
plt.xlabel("Sources $s_1$")
plt.ylabel("Sources $x_2$")
plt.title("Scatter plot of the sources")
fig = plt.gcf()
plt.axis([-1.5*vrange, 1.5*vrange, -1.5*vrange, 1.5*vrange])
fig.set_size_inches(12, 6)
The black arrows display the "mixing directions", in other words the columns of the mixing matrix $\bf A$. When the sources are correctly estimated, these arrows should be aligned with the ordinate and abscissa axes.
BSS is an ill-posed inverse problem, which means that there exists an infinite number of solution factors $({\bf A},{\bf S})$ that can verify the relationship ${\bf X} = {\bf AS}$. Standard approaches try to exploit additional information to limit the number of admissible solutions.
From a very general point of of view, blind source separation has mainly to do with trying to distinguish between the contributions of each individual source. How can we state (and therefore measure) that one source is different from an other ?
Since the beginning of the 90s, the mainstream methods have focused on exploiting statistical properties to measure a difference or diversity between the sources.
I.2 - A first statistical approach: second-order statistics¶
To that respect, the easiest way a statistician can tell that one signal is "different" from an other one is by measuring how thse signal correlates. In other words, in principle, one can easily exploit second-order statistics to mutually characterize the sources: "different" sources are decorrelated signals.
In the next, we will assume that the sources are indeed decorrelated. Moreover, we assume that they are statistical realizations of some stochastic process with mean $0$ and variance $1$. Stating that the variance of the sources is $1$ is not a limitation since the mixture model is invariant according to permutations and to scaling with strictly non-zero factors. This can be formalized as follows:
\begin{array} $\forall i; \; \; \mathbb{E}\{s_i\}&= &0 \\ \forall i; \; \; \mathbb{E}\{s_i^2\}&= &1 \\ \forall i \neq j; \; \; \mathbb{E}\{s_i s_j\}&= &0 \end{array}Looking for sources that are decorrelated can be recast as finding some matrix $\bf A$ such that:
$$ \mathbb{E}\{{\bf X} {\bf X}^T\} = {\bf A} \mathbb{E}\{{\bf S} {\bf S}^T\}{\bf A}^T $$where $\mathbb{E}\{{\bf S} {\bf S}^T\}$ is a diagonal matrix. This therefore amounts to a very basic matrix diagonalization problem, which is usually known as Principal Component Analysis (PCA).
Let's try PCA on the synthetic data.
A_pca,S_pca = pb.Perform_PCA(X,n)
S_pca = np.dot(A_pca.T,X)
S_pca = np.dot(np.diag(1./np.sqrt(np.sum(S_pca*S_pca,axis=1))),S_pca)
vrange = np.max(np.sqrt(np.sum(S_pca*S_pca,axis=0)))
G = np.dot(A_pca.T,A)
G = np.dot(np.diag(1./np.sqrt(np.sum(G*G,axis=1))),G)
G = G.T
plt.title('Sources with PCA')
plt.arrow( 0, 0,vrange*G[0,0], vrange*G[0,1], fc="k", ec="k",linewidth=2,head_width=0.1*vrange, head_length=0.1*vrange)
plt.arrow( 0, 0,vrange*G[1,0], vrange*G[1,1], fc="k", ec="k",linewidth=2,head_width=0.1*vrange, head_length=0.1*vrange)
plt.plot(S_pca[0,:],S_pca[1,:],'co',zorder=0)
plt.axis([-1.5*vrange, 1.5*vrange, -1.5*vrange, 1.5*vrange])
fig = plt.gcf()
fig.set_size_inches(20, 12)
Does it work correctly ? To give a sketch of the answer, let's consider some non-trivial orthogonal matrix $\bf O$ ( i.e. this is basically a rotation matrix that is neither the identity matrix nor a permuation matrix) so that we can define ${\bf A}^\prime = {\bf AQ}$ and ${\bf S}^\prime = {\bf Q}^T {\bf S}$. Then, one can verify that:
These new factors provide a correct reconstruction of the data ${\bf X} = {\bf A}^{\prime} {\bf S}^\prime$.
Additionally, the covariance matrix of the new sources is also the identity matrix.
A as conclusion, looking for decorrelated sources only allows to estimate the sources and the mixing matrix up to an orthogonal or rotation matrix. A complete solution to blind source separation requires a finer statistical characterization of the sources, which allows leveraging the so-called "rotational" degeneracy.
I.3 - From decorrelation to statistical independence¶
Before going any further, let's try to understand how we can discrimimate between components. For that purpose, let's have a look at the scatter plots of mixtures obtained from sources with various statistical distributions: Gaussian, uniform and sparse, after PCA has been applied.
reload(pb)
t = 15000
Xg,Ag,Sg = pb.MakeMixture(SType=1,n=n,m=m,t=t,noise_level = noise_level)
Upg,Xg = pb.Perform_PCA(Xg,n)
Xu,Au,Su = pb.MakeMixture(SType=2,n=n,m=m,t=t,noise_level = noise_level)
Upg,Xu = pb.Perform_PCA(Xu,n)
Xs,As,Ss = pb.MakeMixture(SType=3,n=n,m=m,t=t,noise_level = noise_level)
Upg,Xs = pb.Perform_PCA(Xs,n)
print(np.linalg.norm(Xs[1,:]))
plt.subplot(131)
plt.plot(Xg[0,:],Xg[1,:],'.',color='crimson')
plt.xlabel("Source $x_1$")
plt.ylabel("Source $x_2$")
plt.title("Scatter plot of Gaussian mixtures")
plt.axis([-5/np.sqrt(t), 5/np.sqrt(t), -5/np.sqrt(t), 5/np.sqrt(t)])
plt.subplot(132)
plt.plot(Xu[0,:],Xu[1,:],'.',color='steelblue')
plt.xlabel("Source $x_1$")
plt.ylabel("Source $x_2$")
plt.title("Scatter plot of uniform mixtures")
plt.axis([-5/np.sqrt(t), 5/np.sqrt(t), -5/np.sqrt(t), 5/np.sqrt(t)])
plt.subplot(133)
plt.plot(Xs[0,:],Xs[1,:],'.',color='olive')
plt.xlabel("Source $x_1$")
plt.ylabel("Source $x_2$")
plt.title("Scatter plot of sparse mixtures")
fig = plt.gcf()
plt.axis([-15/np.sqrt(t), 15/np.sqrt(t), -15/np.sqrt(t), 15/np.sqrt(t)])
fig.set_size_inches(15, 5)
To better understand how we can go beyond decorrelation, let's consider that these sources are mixed with orthogonal mixing matrices, in other words rotation matrices $\mathcal{R}({\Theta})$, and let's have a look at the correlation between the sources as the angle $\Theta$ varies.
def RotationMatrix(Theta):
ThetaR = Theta/180.*np.pi
M = [[np.cos(ThetaR),np.sin(ThetaR)],[-np.sin(ThetaR),np.cos(ThetaR)]]
return M
vTheta = np.linspace(-180,180,1000)
sval = np.zeros((3,1000))
cval = np.zeros((3,1000))
for r in range(0,1000):
M = RotationMatrix(vTheta[r])
Xr = np.dot(M,Xg)
cval[0,r] = np.sum(Xr[0,:]*Xr[1,:])
T = np.dot(Xr**3,Xr.T)
sval[0,r] = np.sum(T)
Xr = np.dot(M,Xu)
cval[1,r] = np.sum(Xr[0,:]*Xr[1,:])
T = np.dot(Xr**3,Xr.T)
sval[1,r] = np.sum(T)
Xr = np.dot(M,Xs)
cval[2,r] = np.sum(Xr[0,:]*Xr[1,:])
T = np.dot(Xr**3,Xr.T)
sval[2,r] = np.sum(T)
sval[0,:]= sval[0,:] - np.mean(sval[0,:])
sval[1,:]= sval[1,:] - np.mean(sval[1,:])
sval[2,:]= sval[2,:] - np.mean(sval[2,:])
plt.subplot(121)
plt.plot(vTheta,cval.T,'-',linewidth=2)
plt.xlabel("Rotation angle")
plt.ylabel("Correlation")
plt.title("Correlation of the sources")
plt.legend(['Gaussian','Uniform','Sparse'])
plt.axis([-190, 190, -0.5, 0.5])
plt.subplot(122)
plt.plot(vTheta,sval.T,'-',linewidth=2)
plt.xlabel("Rotation angle")
plt.ylabel("Total skewness")
plt.title("Total skewness of the sources")
plt.legend(['Gaussian','Uniform','Sparse'])
plt.axis([-190, 190, -0.002, 0.002])
fig = plt.gcf()
fig.set_size_inches(20, 10)
As expected, and as shown in the figure on the left, the sources remain decorrelated when mixed with a rotation matrix. However, as shown on the right-hand side pictures, the skewness (statistics of order 3) remains constant for the Gaussian source but varies when the rotation angle varies. On other words, if decorrelation is rotation-invariant, higher-order statistics are not. Therefore, solving BSS problems requires a statistical measure of diversity between the sources that goes beyond second-order statistics ( i.e. 2nd-order statistics provides a Gaussian approximation of the sources' distribution).
In the framework of Independent Component Analysis (ICA), this is done by enforcing the statistical independence between the sought-after sources. If $f_{s_i}$ denotes the probability density function ( p.d.f. ) of source $s_i$ and $f_{\bf S}$ defines the joint p.d.f of the sources, statistical independence is defined by the following relation: $$ f_{\bf S}(s_1,\cdots,s_n) = \prod_{i=1}^n f_{s_i}(s_i), $$ which means that the joint p.d.f of the sources is the product of each marginal p.d.f . In contrast to decorrelation, statistical independence is much more demanding to be numerically enforced. In practice, ICA algorithms are based on indirect measures of this property, whether is comes from statistical approximations/characterization of independence based on higher-order statistics or approaches based on information theory.
Eventually, all these methodes leads to enforcing or emphasizing on the non-Gaussianity of the sources. Supporting this approach, the Darmois-Linnik theorem states that if at most one source is Gaussian and if one finds an unmixing scheme: $$ \hat{\bf S} = {\bf B}{\bf X} $$ so that the sources $\hat{\bf S}$ are statistically independent then the estimated sources are the right sources up to a permutation and a scaling factor.
In this statistical setting, one the fastest/user-friendly generic ICA algorithm is the FastICA algorithm [Hyvarinen99], which enforces the statistical independence of the sources by maximizing their non-Gaussianity. In the following experiments, we will evaluate the performances of the FastICA algorithm.
A_fpica,S_fpica = pb.Perform_FastICA(X,n)
S_fpica = np.dot(np.dot(lng.inv(np.dot(A_fpica.T,A_fpica)),A_fpica.T),X)
S_fpica = np.dot(np.diag(1./np.sqrt(np.sum(S_fpica*S_fpica,axis=1))),S_fpica)
vrange = np.max(np.sqrt(np.sum(S_fpica*S_fpica,axis=0)))
G = np.dot(np.dot(lng.inv(np.dot(A_fpica.T,A_fpica)),A_fpica.T),A)
G = np.dot(np.diag(1./np.sqrt(np.sum(G*G,axis=1))),G)
plt.title('Sources - FastICA')
plt.arrow( 0, 0,vrange*G[0,0], vrange*G[0,1], fc="k", ec="k",linewidth=2,head_width=0.1*vrange, head_length=0.1*vrange)
plt.arrow( 0, 0,vrange*G[1,0], vrange*G[1,1], fc="k", ec="k",linewidth=2,head_width=0.1*vrange, head_length=0.1*vrange)
plt.plot(S_fpica[0,:],S_fpica[1,:],'o',color='khaki',zorder=0)
plt.axis([-1.5*vrange, 1.5*vrange, -1.5*vrange, 1.5*vrange])
fig = plt.gcf()
fig.set_size_inches(12, 12)
Sparse BSS¶
ICA is a very attractive source separation framework with strong theoretical support. However, it suffers from several drawbacks. First, it optimally works if the statistical distribution of the sought-after is known with a certain accuracy, which is barely the case in practice. Second, it generally assumes that the entries/pixels of the sources are independently and identically distributed. However, the precise modelling of natural signal/images generally requires accounting for some correlations between theirs samples as well as non-stationary statistics, which is hard to achieve with standard ICA methods. Based on a radically different approach, more recent BSS methods are grounded upon the sparse modelling of the sources. In a nutshell, the sources are assumed to admit a sparse distribution is some basis/dictionary $\bf \Phi$ ( see the tutorial on sparsity and sparse signal processing ). In the context of BSS, the sparse modelling of the sources is highly beneficial:
- it allows the use simple sparse models to describe potentially highly complex signals/images
- the sparse coding of the information content of the sources allows to better discriminate between different sources, which enhances the separation procedure.
- it provides robustness with respect to non-sparse contaminations such as Gaussian noise
In this framework, the GMCA (Generalized Morphological Component Analysis - [Bobin07]) algorithm aims at solving BSS problems by solving an optimization problem of the form:
$$ \min_{{\bf A},{\bf S}} \sum_i \lambda_i \|s_i {\bf \Phi}^T\|_{\ell_1} + \frac{1}{2}\left \|{\bf X} - {\bf AS} \right\|_F^2, $$where the first (regularization) term enforces sparse sources and the second term measure a discrepancy between the data and the mixture model. The hyperparameters of the model are generalized chosen based on the noise level, which makes it an almost automatic algorithm, which does not require any parameter tuning.
A_gmca,S_gmca,W_gmca = pb.Perform_GMCA(X,n,mints=1,maxts=15)
S_gmca = np.dot(W_gmca,X)
vrange = np.max(np.sqrt(np.sum(S_gmca*S_gmca,axis=0)))
G = np.dot(W_gmca,A)
G = G.T
G = np.dot(np.diag(1./np.sqrt(np.sum(G*G,axis=1))),G)
plt.title('Sources - GMCA')
plt.arrow( 0, 0,vrange*G[0,0], vrange*G[1,0], fc="k", ec="k",linewidth=2,head_width=0.1*vrange, head_length=0.1*vrange)
plt.arrow( 0, 0,vrange*G[0,1], vrange*G[1,1], fc="k", ec="k",linewidth=2,head_width=0.1*vrange, head_length=0.1*vrange)
plt.plot(S_gmca[0,:],S_gmca[1,:],'o',color='steelblue',zorder=0)
plt.axis([-1.5*vrange, 1.5*vrange, -1.5*vrange, 1.5*vrange])
fig = plt.gcf()
fig.set_size_inches(12, 8)
What about noise ?¶
Most statistics-based methods do not explicitly account for the presence of additive noise. However, physical observations are always corrupted with noise. As well, it turns out that the presence of noise can dramatically hamper the performances of blind source separation methods. In the next, we will asumme that the noise $\bf N$ cannot be neglected. For simplicity, we assume that the entries of the noise matrix $\bf N$ are generated independently and identically according to a Gaussian distribution with mean $0$ and some variance $\sigma^2$. Let's start by generating a mixture of sparse sources (2 sources and 2 observations), each one is made of $1024$ samples. The variance $\sigma^2$ is set to that the signal-to-noise ratio is equal to $30$ dB.
n = 2
m = 2
t = 1024
nl = 15 # Signal-To-Noise ratio
SType = 3
X,A,S = pb.MakeMixture(SType=SType,n=n,m=m,t=t,noise_level = nl)
And let's have a look at the scatter of the mixtures when noise is no more negligible.
vrange = 0.5*np.max(np.sqrt(np.sum(X*X,axis=0)))
plt.plot(X[0,:],X[1,:],'.',color='crimson')
plt.xlabel("Data $x_1$")
plt.ylabel("Data $x_2$")
plt.arrow( 0, 0,vrange*A[0,0], vrange*A[1,0], fc="k", ec="k",linewidth=4,head_width=0.1*vrange, head_length=0.1*vrange)
plt.arrow( 0, 0,vrange*A[0,1], vrange*A[1,1], fc="k", ec="k",linewidth=4,head_width=0.1*vrange, head_length=0.1*vrange)
plt.title("Scatter plot of the noisy data")
plt.axis([-2.5*vrange, 2.5*vrange, -2.5*vrange, 2.5*vrange])
fig = plt.gcf()
fig.set_size_inches(20, 12)
Sn = np.dot(np.linalg.inv(A),X)
vrange = 0.5*np.max(np.sqrt(np.sum(Sn*Sn,axis=0)))
plt.subplot(121)
plt.plot(S[0,:],S[1,:],'.',color='crimson')
plt.xlabel("Data $s_1$")
plt.ylabel("Data $s_2$")
plt.title("Scatter plot of the noiseless sources")
plt.subplot(122)
plt.plot(Sn[0,:],Sn[1,:],'.',color='steelblue')
plt.xlabel("Data $s_1$")
plt.ylabel("Data $s_2$")
plt.title("Scatter plot of the noisy sources")
plt.axis([-1.5*vrange, 1.5*vrange, -1.5*vrange, 1.5*vrange])
fig = plt.gcf()
fig.set_size_inches(24, 12)
plt.show
To better measure the impact of noise on the performance of source separation methods, we define comparison criteria:
The mixing matrix criterion: $$ c_A = \frac{1}{n}\left \|{\bf P} \hat{\bf A}^+ {\bf A}\right\|_1, $$ where the matrix $\bf P$ corrects for scaling and permutation. $\hat{\bf A}$ is the estimated mixing matrix and $\hat{\bf A}^+$ is its pseudo-inverse: $\hat{\bf A}^+ = {({\hat{\bf A}}^T\hat{\bf A})}^{-1}{\hat{\bf A}}^T$.
The Signal-To-Distortion Ratio (SDR) measures a discrepancy between the estimated sources and the right sources. It accounts for noise residuals, interferences between badly estimated sources and biais (artifact) of the estimated sources.
The Signal-To-Distortion Ratio (SDR) provides a measure of the interference between the estimated sources.
nc = [40,35,30,25,20,15,10] # Testing for various values of the SNR (in dB)
ntests = 25 # nb of Monte-Carlo simulations
m_fica = np.zeros((len(nc),ntests))
m_gmca = np.zeros((len(nc),ntests))
sdr_fica = np.zeros((len(nc),ntests))
sdr_gmca = np.zeros((len(nc),ntests))
sir_fica = np.zeros((len(nc),ntests))
sir_gmca = np.zeros((len(nc),ntests))
for nt in range(0,ntests):
for r in range(0,len(nc)):
Xo,Ao,So = bu.GenerateMixture(n=n,t=t,m=m,SType=3,noise_level = nc[r])
A_fica,S_fica = pb.Perform_FastICA(Xo,n)
noise = np.dot(np.linalg.inv(A_fica),Xo - np.dot(Ao,So))
SDR,SAR,SIR,SNR = bu.Evaluation_Criteria(Ao,So,A_fica,S_fica)
sdr_fica[r,nt] = SDR
sir_fica[r,nt] = SIR
m_fica[r,nt] = bu.EvalCriterion(Ao,So,A_fica,S_fica)
A_gmca,S_gmca,PiA = bu.GMCA(Xo,n,mints=3,nmax=100,UseP=1,verb=0,percmax=0.5)
S_gmca = np.dot(PiA,np.dot(Ao,So))
noise = np.dot(PiA,Xo - np.dot(Ao,So))
SDR,SAR,SIR,SNR = bu.Evaluation_Criteria(Ao,So,A_gmca,S_gmca)
sdr_gmca[r,nt] = SDR
sir_gmca[r,nt] = SIR
m_gmca[r,nt] = bu.EvalCriterion(Ao,So,A_gmca,S_gmca)
plt.subplot(131)
plt.plot(nc,np.median(m_fica,axis=1),'-d',color='crimson',linewidth=4)
plt.plot(nc,np.median(m_gmca,axis=1),'-o',color='steelblue',linewidth=4)
plt.xlabel("Signal-To-Noise ratio")
plt.ylabel("Mixing matrix criterion")
plt.title("Mixing matrix criterion as a function of the SNR")
plt.legend(["FastICA","Sparse BSS"],bbox_to_anchor=(0., 1.05, 1., .11), loc=3,
ncol=2, mode="expand", borderaxespad=0.)
plt.subplot(132)
plt.plot(nc,np.median(sdr_fica,axis=1),'-d',color='crimson',linewidth=4)
plt.plot(nc,np.median(sdr_gmca,axis=1),'-o',color='steelblue',linewidth=4)
plt.xlabel("Signal-To-Noise ratio")
plt.ylabel("Signal-To-Distortion ratio")
plt.title("SDR as a function of the SNR")
plt.legend(["FastICA","Sparse BSS"],bbox_to_anchor=(0., 1.05, 1., .11), loc=3,
ncol=2, mode="expand", borderaxespad=0.)
plt.subplot(133)
plt.plot(nc,np.median(sir_fica,axis=1),'-d',color='crimson',linewidth=4)
plt.plot(nc,np.median(sir_gmca,axis=1),'-o',color='steelblue',linewidth=4)
plt.xlabel("Signal-To-Noise ratio")
plt.ylabel("Signal-To-Interference ratio")
plt.title("SIR as a function of the SNR")
plt.legend(["FastICA","Sparse BSS"],bbox_to_anchor=(0., 1.05, 1., .11), loc=3,
ncol=2, mode="expand", borderaxespad=0.)
fig = plt.gcf()
fig.set_size_inches(20, 10)
Component separation in a more realistic setting¶
In this section, we will evaluate the performances of these BSS algorihm in a more realistic setting with sources are good approximations of the kind of signals that can be found in spectroscopy. Furthermore, such signals have a sparse distribution in the wavelet domain, which will be exploited to solve BSS problems.
m = 4 # -- defines the number of observations
reload(bu)
# -- Load the sources
X,X0,A0,S0,N,sigma_noise,kern = bu.Make_Experiment_Coherent(n_s=4,n_obs=4,t_samp=1024,w=4,noise_level=120,dynamic=0,sigma1=0.0001,p1=0.001,ptot=0.005)
S0 = abs(S0)
# -- Generate the mixing matrix
A0 = np.random.randn(m,4)
# -- Generate the data matrix
X0 = np.dot(A0,S0)
plt.subplot(121)
plt.title('Spectroscopy - the sources')
plt.plot(S0.T,'-',linewidth=4)
plt.subplot(122)
plt.title('Spectroscopy - the mixtures')
plt.plot(X0.T,'-',linewidth=4)
fig = plt.gcf()
fig.set_size_inches(20, 10)
These sources admit a sparse representation in the wavelet domain. For that purpose, we apply an undecimated wavelet transform to each observation.
J = 3 # -- Number of wavelet scales
reload(bu)
Xw = np.zeros((m,J*1024))
for r in range(0,m):
c,w = bu.Starlet_Forward_1D(X0[r,:],J=J)
Xw[r,:] = w.reshape((1,J*1024))
plt.subplot(121)
plt.title('Scatter plot of 2 mixtures in the sample domain')
plt.plot(X0[0,:],X0[1,:],'o',color=vcol[0],linewidth=4)
vrange = 0.5*np.max(np.sqrt(np.sum(X0[[0,1],:]*X0[[0,1],:],axis=0)))
epsi = 0.02
plt.arrow( 0, 0,vrange*A0[0,0], vrange*A0[1,0], fc="k", ec="k",linewidth=4,head_width=epsi*vrange, head_length=epsi*vrange)
plt.arrow( 0, 0,vrange*A0[0,1], vrange*A0[1,1], fc="k", ec="k",linewidth=4,head_width=epsi*vrange, head_length=epsi*vrange)
plt.arrow( 0, 0,vrange*A0[0,2], vrange*A0[1,2], fc="k", ec="k",linewidth=4,head_width=epsi*vrange, head_length=epsi*vrange)
plt.arrow( 0, 0,vrange*A0[0,3], vrange*A0[1,3], fc="k", ec="k",linewidth=4,head_width=epsi*vrange, head_length=epsi*vrange)
plt.subplot(122)
plt.title('Scatter plot of 2 mixtures in the wavelet domain')
plt.plot(Xw[0,:],Xw[1,:],'o',color=vcol[1],linewidth=4)
fig = plt.gcf()
fig.set_size_inches(20, 10)
We then apply the FastICA algorithm and the GMCA algorithm to these data:
# -- Apply FastICA in the sample domain
A_fpica,S_fpica = pb.Perform_FastICA(X0,4)
S_fpica = np.dot(np.dot(lng.inv(np.dot(A_fpica.T,A_fpica)),A_fpica.T),X0)
S_fpica = np.dot(np.diag(1./np.sqrt(np.sum(S_fpica*S_fpica,axis=1))),S_fpica)
# -- Apply FastICA and GMCA in the wavelet domain
A_fpica_wt,S_fpica_wt = pb.Perform_FastICA(Xw,4)
S_fpica_wt = np.dot(np.dot(lng.inv(np.dot(A_fpica_wt.T,A_fpica_wt)),A_fpica_wt.T),X0)
S_fpica_wt = np.dot(np.diag(1./np.sqrt(np.sum(S_fpica_wt*S_fpica_wt,axis=1))),S_fpica_wt)
A_g,S_g,W_g= bu.GMCA(Xw,4,mints=0,nmax=100,UseP=1,verb=0,percmax=0.5)
S_g = np.dot(lng.inv(A_g),X0)
plt.figure(0)
plt.subplot(131)
plt.title('Scatter plot - FastICA sources in the sample domain')
plt.plot(S_fpica[0,:],S_fpica[1,:],'o',color=vcol[0],linewidth=4)
plt.subplot(132)
plt.title('Scatter plot - FastICA sources in the wavelet domain')
plt.plot(S_fpica_wt[0,:],S_fpica_wt[1,:],'o',color=vcol[1],linewidth=4)
plt.subplot(133)
plt.title('Scatter plot - GMCA sources in the wavelet domain')
plt.plot(S_g[0,:],S_g[1,:],'o',color=vcol[2],linewidth=4)
fig = plt.gcf()
fig.set_size_inches(20, 10)
plt.figure(1)
plt.subplot(131)
plt.title('FastICA sources in the sample domain')
plt.plot(S_fpica.T,linewidth=4)
plt.subplot(132)
plt.title('FastICA sources in the wavelet domain')
plt.plot(S_fpica_wt.T,linewidth=4)
plt.subplot(133)
plt.title('GMCA sources in the wavelet domain')
plt.plot(S_g.T,linewidth=4)
fig = plt.gcf()
fig.set_size_inches(20, 10)
vn = [1,2,3,4,5] # Testing for various values of the peaks' width
ntests = 25 # nb of Monte-Carlo simulations
m_fica = np.zeros((len(vn),ntests))
m_gmca = np.zeros((len(vn),ntests))
sdr_fica = np.zeros((len(vn),ntests))
sdr_gmca = np.zeros((len(vn),ntests))
sir_fica = np.zeros((len(vn),ntests))
sir_gmca = np.zeros((len(vn),ntests))
J = 4
for nt in range(0,ntests):
for r in range(0,len(vn)):
X,Xo,Ao,So,N,sigma_noise,kern = bu.Make_Experiment_Coherent(n_s=5,n_obs=5,t_samp=1024,w=vn[r],noise_level=120,dynamic=0,sigma1=0.0001,p1=0.001,ptot=0.005)
Xw = np.zeros((5,J*1024))
for q in range(0,5):
c,wi = bu.Starlet_Forward_1D(X[q,:],J=J)
Xw[q,:] = wi.reshape((1,J*1024))
A_fica,S_fica = pb.Perform_FastICA(Xw,5)
S_fica = np.dot(np.linalg.inv(A_fica),X)
SDR,SAR,SIR,SNR = bu.Evaluation_Criteria(Ao,So,A_fica,S_fica)
sdr_fica[r,nt] = SDR
sir_fica[r,nt] = SIR
m_fica[r,nt] = bu.EvalCriterion(Ao,So,A_fica,S_fica)
A_gmca,S_gmca,PiA = bu.GMCA(Xw,5,mints=0,nmax=100,UseP=1,verb=0,percmax=0.2)
S_gmca = np.dot(PiA,X)
SDR,SAR,SIR,SNR = bu.Evaluation_Criteria(Ao,So,A_gmca,S_gmca)
sdr_gmca[r,nt] = SDR
sir_gmca[r,nt] = SIR
m_gmca[r,nt] = bu.EvalCriterion(Ao,So,A_gmca,S_gmca)
plt.semilogy(vn,np.median(m_fica,axis=1),'-d',color='crimson',linewidth=4)
plt.semilogy(vn,np.median(m_gmca,axis=1),'-o',color='steelblue',linewidth=4)
plt.xlabel("Signal-To-Noise ratio")
plt.ylabel("Mixing matrix criterion")
plt.title("Mixing matrix criterion as a function of the peaks-width")
plt.legend(["FastICA","Sparse BSS"],bbox_to_anchor=(0., 1.05, 1., .11), loc=3,
ncol=2, mode="expand", borderaxespad=0.)
fig = plt.gcf()
fig.set_size_inches(20, 10)
In the next section, we go on with more complex images. We apply the exact same algorithms to 2D images. We first load some pictures and generate some random mixture.
# --- Load some images
hdulist = pyf.open('./Pictures/AstroImages.fits') # You should specify the path where the data are stored
im = hdulist[0].data
hdulist.close()
J = 3
m = 3
plt.subplot(131)
plt.imshow(im[0,:,:,])
plt.subplot(132)
plt.imshow(im[1,:,:,])
plt.subplot(133)
plt.imshow(im[2,:,:,])
fig = plt.gcf()
fig.set_size_inches(20, 10)
We then apply FastICA and GMCA to these data when the noise level varies.
Np = 64
J = 2
Sim = np.zeros((Np,Np,3))
for r in range(0,3):
Sim[:,:,r] = im[r,0:Np,0:Np]
A0 = np.random.randn(m,3)
U,P,V = np.linalg.svd(A0)
P = np.linspace(0.1,1,3)
A0 = np.dot(U,np.dot(np.diag(P),V.T))
X0 = bu.ApplyAS(A0,Sim)
Xo = np.zeros((3,Np**2))
for q in range(0,3):
Xo[q,:] = X0[:,:,q].reshape(1,Np*Np)
So = np.zeros((3,Np**2))
for q in range(0,3):
So[q,:] = Sim[:,:,q].reshape(1,Np*Np)
isometric = True
wave = RedWave(X0, (0,1), "Coiflet-1",J, True)
# Application of the 2D wavelet transform
Xw = wave.forward(X0)
for scale in range(0,J):
Xw[2*Np*scale:Np + 2*Np*scale,0:Np,:]=0
Xc = np.zeros((3,4*J*Np*Np))
for q in range(0,3):
Xc[q,:] = Xw[:,:,q].reshape(1,4*J*Np*Np)
IndNZ = np.where(abs(Xc[0,:]) > 1e-12)[0]
Sc = np.dot(np.linalg.inv(A0),Xc)
plt.subplot(121)
plt.title('Scatter plot of two mixtures in the wavelet domain')
plt.plot(Xc[0,:],Xc[1,:],'o',color=vcol[4],linewidth=4)
plt.subplot(122)
plt.title('Scatter plot of two sources in the wavelet domain')
plt.plot(Sc[0,:],Sc[1,:],'o',color=vcol[5],linewidth=4)
fig = plt.gcf()
fig.set_size_inches(20, 10)
noise_level = [50,40,30,20]
Ntests = 5
fica_sdr = np.zeros((4,Ntests))
fg_sdr = np.zeros((4,Ntests))
fica_sir = np.zeros((4,Ntests))
fg_sir = np.zeros((4,Ntests))
for nt in range(0,Ntests):
for r in range(0,len(noise_level)):
N0 = np.random.randn(m,Np*Np)
rho = np.power(10,-noise_level[r]/20)*np.linalg.norm(X0)/np.linalg.norm(N0)
X = Xo + rho*N0
A_fpica,S_fpica = pb.Perform_FastICA(X,3)
S_fpica = np.dot(np.linalg.inv(A_fpica),X)
SDR,SAR,SIR,SNR = bu.Evaluation_Criteria(A0,So,A_fpica,S_fpica)
fica_sdr[r,nt] = SDR
fica_sir[r,nt] = SIR
#--- Transform the data in the wavelet domain
Xr = np.zeros((Np,Np,3))
for q in range(0,3):
Xr[:,:,q] = X[q,:].reshape(Np,Np)
Xw = wave.forward(Xr)
for scale in range(0,J):
Xw[2*Np*scale:Np + 2*Np*scale,0:Np,:]=0
Xc = np.zeros((3,4*J*Np*Np))
for q in range(0,3):
Xc[q,:] = Xw[:,:,q].reshape(1,4*J*Np*Np)
A_gmca,S_gmca,W_g= bu.GMCA(Xc[:,IndNZ],3,mints=2,nmax=500,UseP=0,verb=0,percmax=0.2)
S_gmca = np.dot(np.linalg.inv(A_gmca),X)
SDR,SAR,SIR,SNR = bu.Evaluation_Criteria(A0,So,A_gmca,S_gmca)
fg_sdr[r,nt] = SDR
fg_sir[r,nt] = SIR
fg_sdr = np.mean(fg_sdr,axis=1)
fg_sir = np.mean(fg_sir,axis=1)
fica_sdr = np.mean(fica_sdr,axis=1)
fica_sir = np.mean(fica_sir,axis=1)
plt.subplot(121)
plt.plot(noise_level,fg_sdr,color='crimson',linewidth=4)
plt.plot(noise_level,fica_sdr,color='steelblue',linewidth=4)
plt.xlabel('Signal to Noise Ratio')
plt.ylabel('Signal to Distortion Ratio')
plt.legend(['GMCA','FastICA'])
plt.subplot(122)
plt.plot(noise_level,fg_sir,color='crimson',linewidth=4)
plt.plot(noise_level,fica_sir,color='steelblue',linewidth=4)
plt.xlabel('Signal to Noise Ratio')
plt.ylabel('Signal to Interference Ratio')
plt.legend(['GMCA','FastICA'])
fig = plt.gcf()
fig.set_size_inches(20, 10)
noise_level = 20 # Noise level in dB
Np = 256
J = 3
Sim = np.zeros((Np,Np,3))
for r in range(0,3):
Sim[:,:,r] = im[r,0:Np,0:Np]
X0 = bu.ApplyAS(A0,Sim)
Xo = np.zeros((3,Np**2))
for q in range(0,3):
Xo[q,:] = X0[:,:,q].reshape(1,Np*Np)
So = np.zeros((3,Np**2))
for q in range(0,3):
So[q,:] = Sim[:,:,q].reshape(1,Np*Np)
isometric = True
wave = RedWave(X0, (0,1), "Coiflet-1",J, True)
N0 = np.random.randn(m,Np*Np)
rho = np.power(10,-noise_level/20)*np.linalg.norm(X0)/np.linalg.norm(N0)
X = Xo + rho*N0
Xr = np.zeros((Np,Np,3))
for q in range(0,3):
Xr[:,:,q] = X[q,:].reshape(Np,Np)
Xw = wave.forward(Xr)
for scale in range(0,J):
Xw[2*Np*scale:Np + 2*Np*scale,0:Np,:]=0
Xc = np.zeros((3,4*J*Np*Np))
for q in range(0,3):
Xc[q,:] = Xw[:,:,q].reshape(1,4*J*Np*Np)
A_fpica,S_fpica = pb.Perform_FastICA(X,3)
S_fpica = np.dot(np.linalg.inv(A_fpica),X)
S_fpica_wt = np.dot(np.linalg.inv(A_fpica),Xc)
A_gmca,S_gmca,W_g= bu.GMCA(Xc[:,IndNZ],3,mints=2,nmax=500,UseP=1,verb=0,percmax=0.2)
S_g_wt = np.dot(np.linalg.inv(A_gmca),Xc)
S_g = np.dot(np.linalg.inv(A_gmca),X)
plt.subplot(231)
plt.title('FastICA - Source #1')
plt.imshow(S_fpica[0,:].reshape(Np,Np))
plt.subplot(232)
plt.title('FastICA - Source #2')
plt.imshow(S_fpica[1,:].reshape(Np,Np))
plt.subplot(233)
plt.title('FastICA - Source #3')
plt.imshow(S_fpica[2,:].reshape(Np,Np))
plt.subplot(234)
plt.title('GMCA - Source #1')
plt.imshow(S_g[0,:].reshape(Np,Np))
plt.subplot(235)
plt.title('GMCA - Source #2')
plt.imshow(S_g[1,:].reshape(Np,Np))
plt.subplot(236)
plt.title('GMCA - Source #3')
plt.imshow(S_g[2,:].reshape(Np,Np))
fig = plt.gcf()
fig.set_size_inches(30, 20)
PART III - Looking for the CMB in the Planck data¶
You can apply various component separation methods to esimate a CMB map from various regions of the sky.
numb = 1 # --- data to be processed
if numb == 1:
print("PLANCK HFI+LFI - galactic center - 1deg resolution")
hdulist = pyf.open('./CosmoData/Planck_PR1_1deg_galcenter.fits')
PR1 = hdulist[0].data
hdulist.close()
elif numb == 2:
print("PLANCK HFI - galactic center - 10arcmin resolution")
hdulist = pyf.open('./CosmoData/HFI_PR1_10arcmin_galcenter.fits')
PR1 = hdulist[0].data
hdulist.close()
elif numb == 3:
print("PLANCK HFI - Taurus molecular cloud - 10arcmin resolution")
hdulist = pyf.open('./CosmoData/HFI_PR1_10arcmin_taurus_MC.fits')
PR1 = hdulist[0].data
hdulist.close()
if numb == 1:
plt.figure(0)
plt.subplot(331)
plt.title('30 GHz')
plt.imshow(PR1[0,:,:],vmin=-5*pb.mad(PR1[0,:,:]),vmax=5*pb.mad(PR1[0,:,:]))
plt.subplot(332)
plt.title('44 GHz')
plt.imshow(PR1[1,:,:],vmin=-5*pb.mad(PR1[1,:,:]),vmax=5*pb.mad(PR1[1,:,:]))
plt.subplot(333)
plt.title('70 GHz')
plt.imshow(PR1[2,:,:],vmin=-5*pb.mad(PR1[2,:,:]),vmax=5*pb.mad(PR1[2,:,:]))
plt.subplot(334)
plt.title('100 GHz')
plt.imshow(PR1[3,:,:],vmin=-5*pb.mad(PR1[3,:,:]),vmax=5*pb.mad(PR1[3,:,:]))
plt.subplot(335)
plt.title('143 GHz')
plt.imshow(PR1[4,:,:],vmin=-5*pb.mad(PR1[4,:,:]),vmax=5*pb.mad(PR1[4,:,:]))
plt.subplot(336)
plt.title('217 GHz')
plt.imshow(PR1[5,:,:],vmin=-5*pb.mad(PR1[5,:,:]),vmax=5*pb.mad(PR1[5,:,:]))
plt.subplot(337)
plt.title('353 GHz')
plt.imshow(PR1[6,:,:],vmin=-5*pb.mad(PR1[6,:,:]),vmax=5*pb.mad(PR1[6,:,:]))
plt.subplot(338)
plt.title('545 GHz')
plt.imshow(PR1[7,:,:],vmin=-5*pb.mad(PR1[7,:,:]),vmax=5*pb.mad(PR1[7,:,:]))
plt.subplot(339)
plt.title('857 GHz')
plt.imshow(PR1[8,:,:],vmin=-5*pb.mad(PR1[8,:,:]),vmax=5*pb.mad(PR1[8,:,:]))
else:
plt.figure(0)
plt.subplot(231)
plt.title('100 GHz')
plt.imshow(PR1[0,:,:],vmin=-5*pb.mad(PR1[0,:,:]),vmax=5*pb.mad(PR1[0,:,:]))
plt.subplot(232)
plt.title('143 GHz')
plt.imshow(PR1[1,:,:],vmin=-5*pb.mad(PR1[1,:,:]),vmax=5*pb.mad(PR1[1,:,:]))
plt.subplot(233)
plt.title('217 GHz')
plt.imshow(PR1[2,:,:],vmin=-5*pb.mad(PR1[2,:,:]),vmax=5*pb.mad(PR1[2,:,:]))
plt.subplot(234)
plt.title('353 GHz')
plt.imshow(PR1[3,:,:],vmin=-5*pb.mad(PR1[3,:,:]),vmax=5*pb.mad(PR1[3,:,:]))
plt.subplot(235)
plt.title('545 GHz')
plt.imshow(PR1[4,:,:],vmin=-5*pb.mad(PR1[4,:,:]),vmax=5*pb.mad(PR1[4,:,:]))
plt.subplot(236)
plt.title('857 GHz')
plt.imshow(PR1[5,:,:],vmin=-5*pb.mad(PR1[5,:,:]),vmax=5*pb.mad(PR1[5,:,:]))
fig = plt.gcf()
fig.set_size_inches(20, 10)
c = np.array([0.97707402,0.95145921,0.88249613,0.77729534,0.60483320,0.33441687,0.077544524,0.0062673491,6.3773948e-05])
X = bu.Cube2Mat(PR1)
nobs,nx,ny = np.shape(PR1)
Internal Linear Combination - ILC¶
ILC is CMB estimation technique that is very popular in astrophysics. In a nutshell, it performs by estimating a single source $s$ with known spectrum $a$ by minimizing its variance: $$ \min_{w; \, w^T a =1} \mbox{Var}\{ w {\bf X}\} $$ where the source $s$ is defined as a linear combination of the input observations with some weights $w$.
if numb == 1:
colcmb = c
else:
colcmb = c[3:9]
s_ilc,w_ilc = pb.Perform_ILC(X,colcmb)
s_ilc = s_ilc - np.mean(s_ilc)
plt.subplot(131)
plt.imshow(np.reshape(s_ilc,(nx,ny)),vmin=-3*pb.mad(s_ilc),vmax=3*pb.mad(s_ilc))
plt.title('CMB with ILC')
fig = plt.gcf()
fig.set_size_inches(20, 10)
FastICA¶
In this section, we apply the FastICA algorithm to the Planck data.
A_fpica,S_fpica = pb.Perform_FastICA(X,nobs)
for r in range(nobs):
rho = lng.norm(colcmb)/lng.norm(A_fpica[:,r])
A_fpica[:,r] = A_fpica[:,r]*rho
S_fpica[r,:] = S_fpica[r,:] / rho
G_fpica = np.dot(lng.inv(A_fpica),np.reshape(colcmb,(len(colcmb),1)))
index_cmb = np.argsort(abs(G_fpica),axis = 0)[::-1]
s_sign = np.sign(G_fpica[index_cmb[0]])
s_fpica = s_sign*S_fpica[index_cmb[0],:]
plt.figure(1)
plt.subplot(132)
plt.imshow(np.reshape(s_fpica,(nx,ny)),vmin=-3*pb.mad(s_fpica),vmax=3*pb.mad(s_fpica))
plt.title('CMB with FastICA')
fig = plt.gcf()
fig.set_size_inches(20, 10)
plt.figure(2)
plt.subplot(131)
plt.imshow(np.reshape(s_fpica,(nx,ny)),vmin=-3*pb.mad(s_fpica),vmax=3*pb.mad(s_fpica))
plt.title('1st FastICA component')
plt.subplot(132)
plt.imshow(np.reshape(S_fpica[index_cmb[1],:],(nx,ny)),vmin=-3*pb.mad(s_fpica),vmax=3*pb.mad(s_fpica))
plt.title('2nd FastICA component')
plt.subplot(133)
plt.imshow(np.reshape(S_fpica[index_cmb[2],:],(nx,ny)),vmin=-3*pb.mad(s_fpica),vmax=3*pb.mad(s_fpica))
plt.title('3rd FastICA component')
fig = plt.gcf()
fig.set_size_inches(20, 10)
GMCA¶
In this section, we apply the GMCA algorithm to the Planck data.
if numb == 1:
hdulist = pyf.open('./CosmoData/WT_Planck_PR1_1deg_galcenter.fits')
wt_PR1 = hdulist[0].data
hdulist.close()
elif numb == 2:
hdulist = pyf.open('./CosmoData/WT_HFI_PR1_10arcmin_galcenter.fits')
wt_PR1 = hdulist[0].data
hdulist.close()
elif numb == 3:
hdulist = pyf.open('./CosmoData/WT_HFI_PR1_10arcmin_taurus_MC.fits')
wt_PR1 = hdulist[0].data
hdulist.close()
fig.set_size_inches(20, 10)
wt_X = pb.Cube2Mat(wt_PR1)
A_g,S_g,W_g= pb.Perform_GMCA_PLANCK(wt_X,nobs,nmax = 100,mints = 1,colcmb=colcmb)
S_g = np.dot(lng.inv(A_g),X)
plt.figure(1)
plt.subplot(133)
plt.imshow(np.reshape(S_g[0,:]-np.median(S_g[0,:]),(nx,ny)),vmin=-3*pb.mad(S_g[0,:]),vmax=3*pb.mad(S_g[0,:]))
plt.title('CMB with GMCA')
fig = plt.gcf()
fig.set_size_inches(20, 10)
Diff = S_g[0,:] - s_ilc
plt.figure(3)
plt.subplot(121)
mx=3*pb.mad(s_ilc)
plt.imshow(np.reshape(s_ilc-np.median(s_ilc),(nx,ny)),vmin=-mx,vmax=mx)
plt.colorbar(shrink=0.5)
plt.title('ILC map')
plt.subplot(122)
plt.imshow(np.reshape(Diff-np.median(Diff),(nx,ny)),vmin=-0.25*mx,vmax=0.25*mx)
plt.colorbar(shrink=0.5)
plt.title('Difference between GMCA and ILC')
fig = plt.gcf()
fig.set_size_inches(20, 10)
We advise you to try these codes on different regions of the sky. You can clearly observe that all these methods yield to roughly acceptable estimates of the CMB but do not provide perfectly clean pictures.
The main limitation here is that the linear mixture model ${\bf X} = {\bf AS} + {\bf N}$ is not perfectly valid with these data. In fact, this model states that the contribution of each astrophysical component is a rank-1 matrix, which is composed of a single electromagnetic spectrum (column of the matrix $\bf A$) and a spatial distribution (i.e. the source). This model is a very good approximation for components like the CMB or the Sunyaev-Zel'Dovich effect but not for galactic foregrounds such as synchroton, free-free and thermal dust emission.
BIBLIOGRAPHY¶
- [Hyvarinen99] "Fast and Robust Fixed-Point Algorithms for Independent Component Analysis" , IEEE Transactions on Neural Networks, 10(3), 1999.
- [Cardoso98] "Blind Signal Separation: Statistical Principles" , Proceedings of the IEEE, 9(10), 1998.
- [Comon10] "Handbook of blind source separation" , Academic Press, 2010.
- [Cichocki07] "Sparse Component Analysis" , Data Mining in Biomedicine, Springer, 7, 2007.
- [Bobin07] "Sparsity and Morphological Diversity in Blind Source Separation " , IEEE Transactions on Image Processing, 16(11), 2007.
- [Bobin13] "Sparse component separation for accurate CMB map estimation" , Astronomy and Astrophysics, 550(A73), 2013.
- [Bobin14] "Joint Planck and WMAP CMB map reconstruction " , Astronomy and Astrophysics, 563(A105), 2014.
- [Bobin15] "Sparsity and adaptivity for the blind separation of partially correlated sources" , IEEE Transactions on Signal Processing, 63(5), 2015.
- [Starck07] "The Undecimated Wavelet Decomposition and its Reconstruction" , IEEE Transactions on Image Processing, 16(2), 2007.
- [Starck10] "Sparse Image and Signal Processing" , Cambridge University Press, 2010.